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The radial Schrodinger equation for a spherically symmetric potential can be regarded as a one 
dimensional classical harmonic oscillator with a time-dependent spring constant. For solving classical 
dynamics problems, symplectic integrators are well known for their excellent conservation properties. 
The class of gradient symplectic algorithms is particularly suited for solving harmonic oscillator 
dynamics. By use of Suzuki's rule for decomposing time-ordered operators, these algorithms can 
be easily applied to the Schrodinger equation. We demonstrate the power of this class of gradient 
algorithms by solving the spectrum of highly singular radial potentials using Killingbeck's method 
of backward Newton-Ralphson iterations. 

I. INTRODUCTION 

Because of its physical importance, an immense literature exists for solving the radial Schrodinger equation, 

d 2 u(r) 



dr 2 

where 



f(r,E)u(r), (1.1) 



f(r,E) = 2V(r)-2E+ £ -^±^-. (1.2) 

This is usually solved by finite difference methods, such as the well known fourth-order Numero\si algorithm, or 
further improved schemes^. Recent devlopments^i^ have resulted in many exponentially fitted algorithms which seek 
to integrate (|l.lf) exactly when / is a constant. As we will see, because / can vary rapidly with V(r), specially in the 
case of singular potentials, these algorithms do not in general perform better than non-fitted algorithms. 

If we relabel the variables r — > t and u — > q, l|l.l|) is just a ID harmonic oscillator with a time-dependent spring 
constant k{t,E) = —f(t,E), 

fg±=- k (t,E) q (t). (1.3) 

The novel aspect of the problem is that k(t,E) can change sign with time and become repulsive beyond the turning- 
point. The dynamics of l|1.3|) corresponds to a Hamiltonian with an explicit time-dependent potential, 

H= l -p 2 + l -k{t,E)q 2 . (1.4) 

Thus any algorithm that can solve the classical time-dependent force problem can be used to solve the radial 
Schrodinger equation. For example, one can use Runge-Kutta type algorithms^. However, for solving classical 
dynamics, symplectic integratorsi&Siifi are algorithms of choice because they conserve all Poincare invariants and 
are deeply rooted in the Poisson formulation of classical mechanics. For oscillatory problems, symplectic algorithms 
are known to conserve energy and reduce phase error much better than Runge-Kutta type algorithmsiiiiSiiiiii^. 
The difficulty here is that, in order to derive an algorithm for solving time-dependent dynamics, one must solve the 
problem of time-ordered exponential. Liu et nlM± have recognized the time-dependent Hamiltonain structure of the 
Schrodinger equation, but were able to solve the time-dependent exponential, and devised a symplectic algorithm, 
only to second order. Kalogiratou, Monovasilis and Simosi£ have proposed a third order symplectic algorithm by 
expanding out the exponential to third order. Such a brute force approach cannot be extended to higher orders. 
A more systematic way of dealing with the time-ordered exponential is via the Magnus expansioni£ii2*2ii, but the 
Magnus expansion requires explicit time integration in addition to evaluating higher order commutators. A more ele- 
gant solution is Suzuki'sSi reinterpretation of the time-order exponential as reviewed in Ref^. By adapting Suzuki's 
rule, any factorized symplectic algorithms can be used to solve problems with an explicit time-dependent potential^, 
including the disguised radial Schrodinger equation <|1.3|) . 

In order to devise efficient algorithms for solving the radial Schrodinger equation i|l.l[l . one must take advantage 
of its harmonic oscillator character l|1.3f) . Most algorithms, even factorized symplectic ones, are general purpose 
algorithms and are not specially adapted for solving the time-dependent harmonic oscillator. However, the recent 
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class of gradient symplectic algorithmsi2ii&i»ii£i22i2ii2», while general, seem tailor-made for solving harmonic type 
dynamics. This is because these algorithms require computing the force gradient in addition to the force. While the 
force gradient is not difficult to compute, it is especially trivial in the case of the harmonic oscillator. This class of 
gradient (or more specifically, forward) integrators, has been shown to be efficient in solving both classicaliSii&iiii^iSiSi 
and quant umi2i22i22i2SiS& dynamical problems. In this work, we will show they are also ideally suited for solving the 
radial Schrodinger equation. 

In the next section, we briefly summarize the Lie-Poisson operator formulation of symplectic integrators and Suzuki's 
rule for factorizing time-ordered exponentials. In Section III, we describe forward, gradient based symplectic algo- 
rithms. In Section IV we review Killingbeck's method^ 7 of eigenvalue-function determination. In Section V, we 
compare results on the Coulomb and other singular radial potentials. In Section VI, we discuss the applicability of 
sixth-order algorithms and draw some conclusions in Section VII. 



II. TIME-DEPENDENT SYMPLECTIC ALGORITHMS 



The Poisson bracket for evolving any dynamical variable W{q,p) can be regarded as an operator equation 

| CD 

with formal solution 

W(t + e) = c e{T+v ^W{t). (2.2) 

For the standard Hamiltonian, 

H(p,q) = ^p 2 + V(q), (2.3) 

the operators T and V are first order differential operators 

T _dH_d__ d_ __d£d__ d_ 

dp dq 9g' dq dp ^ dp' 

The Lie transforms^ e eT and e eV , are then displacement operators which displace q and p forward in time via 

q — > q + sp and p — > p + eF. (2.5) 

Each factorization of e £ ( T+v * 1 into products of e eT , e eV (and exponentials of commutators of T and V) gives rise 
to a symplectic algorithm, which is a sequence of successive displacements l|2.5|l for evolving the system forward in 
time. This is the fundamental Lie-Poisson theory of symplectic integrators which has been studied extensively in the 
literaturo 7 ! 8 ! 9 ! 10 . 

For a time-dependent Hamiltonian 

H{t)= l -p 2 + V{q,t), (2.6) 
the solution is given by the time-order exponential 

W(t + e) = Texp(J [T + V(s)]dsjW(t), (2.7) 
where V(t) is now the explicitly time-dependent operator 

Suzuki proved^i that 

Texp([ t+6 {T + V( S )}ds) =C WM+^ (2 . 9) 
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where D is the forward time derivative operator 

D = M ^ 

such that for any two time-dependent functions F(t) and G(t), 

F(t)e eD G(t) = F{t + e)G(t). (2.11) 

Thus symplectic algorithms for solving explicitly time-dependent problems of the form (|2.6(l can be obtained by 
factorizing the three-operator exponential of (|2.9[1 . Since D commutes with T, one can first group 

f = T + D (2.12) 

and factorize T and V(t) as in the time- independent case. The difference between time-dependent and time- 
independent algorithms resides solely in the use of T in place of T. This makes it extremely easy to analyze and 

devise time-dependent algorithms. Once factorized in terms of c eT = c sT c eD , the operator e eD then shifts the time 
at which all the time-dependent potential to its LEFT must be evaluated. This results in Suzuki's rule for solving 
time-dependent Hamiltonian Ij2.6|) : the time- dependent potential must be evaluated at an incremental time from the 
start of the algorithm equal to the sum of time steps of all the T operators to its right. We will illustrate how this is 
applied in the next section. For more detailed discussions and examples, see Refsi 14 i 22 . 

III. FORWARD FOURTH-ORDER ALGORITHMS 

In order to solve the radial Schrodinger equation l|l.l|l efficiently, one must take advantage of its harmonic oscillator 
character (jl.3[) . This can be done easily for factorized algorithms because their error terms have a well-defined 
analytical structure. Consider the second order factorization, 



-eT eV -eT 

e 2 e e 2 = expe 



(T + V) + -e 2 ([T, [V, T]} - 2[V, [T, V]]) + 0(e 4 ) \ . (3.1) 



This is just a general operator equality which follows from the Baker-Campbell-Hausdorff (BCH) formula. In the 
present context, this equlity tells us that the second-order factorization on the LHS deviates from the exact evolution 
operator exp e ( T+v ^ by error terms which are the double commutators on the RHS. However, for the ordinary harmonic 
oscillator Hamiltonian l|1.4|) with a constant spring constant k — ui 2 , one can easily verify that 

[V,[T,V]] = -2lu 2 V, (3.2) 

[T, [V, T]] = -2lu 2 T. (3.3) 

Thus the error terms can be re-expressed in terms of the original operator T and V and be moved back to the LHS 
to yield, 

e e(I + Jj^ 2 £ 2 )T e e(l-i^ 2 e 2 )y e £(i + ^^ 2 e 2 )T = ^(T+V +0(E i )) _ (3 4 ) 

This means that the LHS is now a fourth-order algorithm for solving the harmonic oscillator. Because of the funda- 
mental identities l|3.2|) and (|3.3I) , all higher order commutators for the harmonic oscillator can be subsummed back to 
T and V yielding the exact factorization 26 

e EC E T e EC M V e EC E T = e e(T+V) ; ( g 5) 

where the "edge" cocfficcint Cg and the "middle" coefficient Cm are given by 

1 - cos(we) sin(we) 
C E = — — - and C M = • (3.6) 

we sin(w£) cue 

The above discussion only depends on the abstract commutator relations (|3.2|) and l|3.3|l . and is independent of the 
specific form of the operator T and V. Thus by interchanging T «-> V , we can also factorize exactly, 

e e{T+V) = e EC E V e EC M T e EC E V _ (37) 
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To solve the time-dependent harmonic oscillator, one has to replace T — > T everywhere. It is easy to verify that for 
any two time-dependent functions W(t) and V(t), 

[V(t),[D,W(t)}}=0. (3.8) 

Hence, the commuatator (|3.2H 

[V, [f, V}} = [V, [T, V}} = 2f(t, E)V (3.9) 

remains proportional to V. However, 

[T, [V, T]] = 2/(t, E)T + 2[T, -V] - (3.10) 

bears no simple relationship to T. This means that one can still retain the commutator [V, [T, V]] — [V, [T, V]] and 
move it back to the LHS, but in order to have a fourth-order algorithm, one must eliminate the commutator [T, [V, T]] 
by more elaborated factorization schemes. This coincides precisely with the way forward symplectic algorithms are 
derivedi 2 *^* 2 ^ For example, the simplest fourth-order forward factorization scheme 12 ' 28 4B for evolving the system 
forward for time e is 

Xj 4) ( £ ) = e a£5 V^*e te? e £ * y *e° £? , 

= e asT e e±V*(t+a's) e bsT e e±V(t+ae) e aeT ^ /g -qs 

where a = |(1 — ^=), b = -^=, a' = a + b = i(l + -i=), and the effective potential operator is given by 

V*(t) = V(t) + 1(2 - V3)s 2 [V(t), [T, V(t)]] 

l + l(2-^)e 2 /(t,i?)]^). (3.12) 

This results in the use of an effective time-dependent force 

F*(t, E) = [1 + 1(2 - V3)e 2 /(t, £0] /(*, E)g, (3.13) 

which is no more difficult to evaluate than the original. Factorization scheme (I3.11|) translates into the following 
fourth-order algorithm for solving the time-dependent harmonic oscillator: 

<?i = Qa + ae p 

Pi = Po + -^£F*(t + ae)q 1 

Q2 = qi+ bepi (3.14) 
P2 = pi + ^sF*(t + a'e)q 2 
q3 = q2 + ae p 2 

The last numbered p and q are the updated values. In the present context, since q is the wave function and p is 
only an ancillary variable, we will be interested only in algorithms that begins and ends with q. These position-type 
algorithms make full use of force evaluations at intermediate time to update the final position. As will be discuss in 
the next Section, this point is important for Killingbeck's method of iterating the last position q(0, E) to zero. 
In general, the commutator 



dFi d „ d 

3~ 



[V, [T, V}} = 2F,p-^ = V,(|F| 2 )-^ (3.15) 
dqj dp l dpi 



produces a force which is the gradient of the square of the original force. For the ID harmonic oscillator, this is simply 
d(fq) 2 /dq = 2f 2 q. By incorporating the force gradient, algorithm 4B (|3.14() is fourth-order with only two evaluations 
of the effective force Ij3.13ll . 

For three force evaluations, one can use algorithm 4G— : 

T^ 4) (e) = e i sT ei eV+a ^ e3u e^ eT ei eV+ ^- 2a ^ e3u ei sT ei eV+a ^ s3u ei eT , (3.16) 
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where U = [V, [T, V]]. One is free to distribute the commutator term symmetrically via a without affecting its 
fourth-order convergence. The three obvious choices are a — 0, 3/8, 1/2. The first and the last case concentrate the 
gradient term at the center and at the two sides respectively. The second case distributes the gradient term in the 
same proportion as the original force so that the same effective force 

F*(t, E) = [l+ ^£ 2 f(t, E)] f(t, E)q, (3.17) 

is evaluated at three different times. This is a direct generalization of algorithm 4B. We shall refer to these three 
variants as 4C, 4C' and 4C". For any specific application, one can fine-tune a to minimize, or even eliminate, the 
algorithm's fourth-order step-size error. We shall refer to this optimized case as 4C op t. Other forward, or just gradient 
based algorithms, can be found in 12 i 14 i 22 i 23 i 24 . 



IV. KILLINGBECK'S BACKWARD ITERATION 



Killingbeck's method^ 7 , for solving the eigenvalue- function pair requires no wave function matching and can be highly 
automated. It consists of two key steps: 1) backward integration to ensure numerical stability, and 2) quadratic energy 
convergence via Newton-Ralphson iterations. One begins with an initial guess of the eigenvalue E 1 - ^ and chooses a 
large time value T (large R in the original problem) to set q(T) = and p(T) = p^, where is an arbitrary, but 
small number. One then iterates the algorithm, such as l|3.14|l . backward in time to t = 0. (In practice, it may be 
simpler to run the algorithm forward in time and change the potential from V(t) to V(T — t).) If E is a correct 
eigenvalue, then it must satisfy the eigencondition 

q(0,E) = 0. (4.1) 
Thus the eigenvalue E is a root of the above equation and can be solved by Newton-Ralphson iterations: 

Killingbeck suggested that the derivative q'(Q, E) = dq(Q, E)/dE, which obeys, 

<fq'{t) 



dt 2 



= f(t,E)q'(t)-2q(t), (4.3) 



can be solved simultaneously with q(t), i.e., differentiating any algorithm, such as l|3.14|l . line by line with respect to 
E. The resulting algorithm can be iterated at the same time to determine both q(0,E) and q'(0,E) simultaneously 
so that 1)4. 2J) can be updated directly. By re-running the algorithm with the updated energy, the procedure can be 
repeated until convergence. The convergence is quadratic in the number of iterations. The converged eigenvalue 
(and eigenfunction) will deviate from the exact value in powers of e" depending on the order of the algorithm used. 
In solving the radial Schrodinger equation, t — (i.e., r — 0) is the absolute boundary and f(t,E) is not defined 
for t < 0. Thus in applying Killingbeck's method, one must not use any algorithm which evaluate the force at an 
intermediate time greater than t + e. 



V. RESULTS FOR SINGULAR POTENTIALS 



One important application of solving the radial Schrodinger equation is in atomic {e.g., density functional) calcu- 
lations, where the dominant interaction is the Coulomb potential 

V{r)=-- (5.1) 
r 

As a prototype test case, we show the convergence of various algorithms in solving for the ground state of (|5.1|l in 
Figure 1. We use T = 26 (R = 26); beyond T — 25, there is no change in the eigenvalue on the order of 1CP 12 . For an 
initial guess of E^ = —0.6, the Killingbeck iteration converges to 12 decimal places in 9 iterations or less. In most 
cases, once a good guess is found, only a few iterations are necessary. 

It is well known that when the Numerov (N) algorithm 1 , is used in Killingbeck's method, the Coubomb ground 
state only converges quadraticallj* 2 ^. While the reason for this is understood-S and a simple remedy is available 31 , 
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most of the self-starting fourth-order algorithms used here suffered no such order reduction. Most can be well-fitted 
by the power law E(e) — E = ce 4 . RKN and RK are the three and four force-evaluation Runge-Kutta-Nystrom 
and Runge-Kutta algorithms^, respectively. FR is the Forest-Ruth^ symplectic algorithm which uses three force- 
evaluations. This is the first fourth-order symplectic algorithm found and is well known for its relative large error. 
M is McLachlan's improved fourth-order algorithmic which uses four force-evaluations. BM is Blanes and Moan's 
latest^ refined fourth-order algorithm which uses six force-evaluations. 

FR, M and BM are examples of conventional symplectic algorithms which have negative intermendiate time-steps. 
As shown in Figure 1, algorithm 4B, which use only two evaluations of the effective force, outperforms all the 
aforementioned algorithms except BM regardless the number of force-evaluation. OMF18, OMF29, and OMF36 are 
Omelyan, Mryglod and Folk's listed^ fourth-order algorithms 18, 29 and 36. These are gradient algorithms, similar 
to 4B and 4C', which use three, four, and five effective force evaluations, respectively. As a is varied from to 0.5, 
the error of the general 4C algorithm changes from negative to positive. At a = 0.49, the error curve resembles that of 
BM. At the optimal value of a = 0.41, the fourth-order error should have been nearly eliminated with the algorithm 
showing sixth-order convergence. The fact that it does not will be discussed in the next section. We fitted all the 
results in Figure 1 via a power-law of the form E(e) — Eq = ce n to verify their order of convergence. All can be 
well-fitted with n = 4 except 4B and 4C opt at a = 0.41. For 4B, n w 3.5 and for 4C opt , n w 4.5. Why algorithm 4B 
should suffer such an order reduction is not understood. It is possible that for 4B, its power-law behavior only sets 
in at smaller e. The case of 4C op t will be discuss in the next section. 

Since algorithm OMF29 uses four effective force evaluations, one can run algorithm 4B twice at the half the 
time step size. Thus one should compare OMF29 with 4B(e/2), or OMF29(2e) with 4B(e). Thus relative to the 
computational effort of4B(e), one should compare 4C'(1.5e), OMF18(1.5e), OMF29(2e), OMF36(2.5e) and BM(3e). 
This comparison is shown in Figure 2. In this equal effort comparison, algorithm 4B's fourth order error is as small, 
if not smaller than all the other gradient algorithm's error. This illustrates the case that efficiency is not necessarily 
enhanced by increasing the number of force evaluations. Also, all gradient algorithms have errors smaller than that 
of BM despite fewer force evaluations. We conclude that in solving the Coulomb ground state, the efficiency of 
algorithms 4B and 4C ' is unsurpassed by any other algorithms except by the tunable 4C algorithm. 

In the first column of Table I, we list the energy obtained by all the algorithms at e = 0.01 weighted by their number 
of force evaluations. Algorithm 4B and 4C ' indeed turn in the best result and are outperformed only by 4C op t at 
a = 0.41. For more accurate results, one can just reduce e. 

As a more stringent test of our algorithms and Killingbeck's method, we next consider the spiked harmonic oscillator 
(SHO) with potential 

V(r) = l(r 2 + ±y (5.2) 

For extensive references and discussion on SHO, see Refs^^SiSLS. Figure 3 shows the convergence of the ground 
state energy for the well studied case2£i2L2& of M = 6 with A = 0.001. For T = 10 and a reasonable initial guess of 
E^ ' — 1.5, only five or less iterations are needed for the energy to converge to 12 decimal places. For such a singular 
potential, the convergent step size has to be much smaller, but surprisingly, only a magnitude smaller. Despite the 
high degree of singularity, nearly all algorithms remained fourth order and none is down graded to lower order. At 
a glance, all gradient based algorithms converge better than non-gradient algorithms. Even BM is no better than 
algorithm 4C. Since 4C and 4C" have errors of opposite sign, one can again vary a to minimize the fourth order 
error. The optimized algorithm at a = 0.22 converges better than all other algorithms regardless of the number of 
force evaluations and can be best fitted by a fifth order power law. 

To compare the efficiency of gradient algorithms, we again normalize each algorithm to the computation effort 
of 4B. In Figure4, we plot the convergence curve of 4B(e), 4C'(1.5e), OMF29(2e) and OMF36 (2.5 e). The solid 
line is the fourth order monomial E(e) — Eq = ce 4 which goes through 4C"s result with c = 16.7. The other three 
algorithms can be fitted with the dotted line with c = 20.0. Thus all gradient algorithms are essentially similar, with 
4C marginally better. Again, algorithms OME29 and OME36, which use four and five effective force evaluations 
with complex numeric coefficients, are not more efficient than the simpler algorithms 4B and 4C ' with analytical 
coefficients. 

At e = 0.001, the weighted result of each algorithm is given in the second column of Table I. All are in excellent 
agreement with the value found in the literatureS^LiS. At this step size, only gradient algorithms are accurate to 10 
or more decimal places. For even greater accuracy, one can simply reduce e. 

The algorithms are equally effective in the case of M = 4 and A = 0.001. This is shown in Figure 5. All algorithms 
showed fourth order convergence, except for 4C op t, which can be better fitted with a fifth order power-law. Their 
energy values at e = 0.001 are listed in the third column of Table I. The optimized 4C algorithm is accurate to 9 
decimal places. Note that once algorithm 4C is optimize for M = 6, it can also be used for M = 4. The change in 
a's value is slight. 
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In the most difficult, "supersingular" case of M — 5/2, with A remained small at 0.001, the power-law behavior 
seems to require e < 10~ 5 . This is shown in Figure 6. (If A were not too small, such as 0.1 or 0.01, the power-law 
behavior would remain observable in the range of e considered.) The energy obtained at e = 0.0002 is listed in the 
fourth column of Table I. The variable 4C algorithm uses a — 0.23 inherited from the M = 4 case. All algorithms are 
less efficient in dealing with this "supersingular" case, but gradient algorithms can still maintain an 8-digit accuracy. 

VI. HIGHER ORDER ALGORITHMS 

In general, if Ta is a left-right symmetric approximation of the short time evolution operator e e ^ T+v \ 

N 

T A =Y\c U£T e v * eV = e £liA , (6.1) 

such that ti — 1 and Vi = 1, then the approximate Hamiltonian operator is of the form 

H A = T + V + e 2 ( e TTV [T 2 V] + e VTV [VTV]) 

+ £ 4 ( e TTTTV [T T 3 V] + e VTTTV [V T 3 V] 

+ e TTVTV [TTVTV] + e VTVTV [VTVTV}) + ..., (6.2) 

where e TTV , e VTTTV etc., are coefficients specific to a particular algorithm and where we have used the condensed 
commutator notation [T 2 U] = [T, [T, V]], etc.. Symmetric factorizations give rise to time-reversible algorithms and 
have only even-order error terms. For a constant k = ui 2 , the fundamental identity 13.311 implies that [TT 3 U] = 
and |VT 3 V"] = 0. This crucial simplification is no longer true in the time-dependent case when T is replaced by 
T. From this perspective, one can understand why the ability to integrate the time-independent harmonic oscillator 
exactly does not help in solving the time-dependent case. Exponentially or sinusoidally fitted algorithms are therefore 
not necessarily more efficient. In the time-dependent case, the problem is fundamentally different because some 
commutators no longer vanish. Note that the commutators 

[VTVTV] = 4:f 2 (t)V (6.3) 

can be moved back to the LHS. However the saving here is marginal since this error term can be easily eliminated by 
incorporating more gradient terms [VTV] in the factorization process. 

It has been showniS that in order to derived a general sixth-order forward algorithm, one must retain both [T^TU] 
and [VT 3 V] in the factorization process. Unfortunately, since [UT 3 U] cannot be evaluated easily, there is currently no 
practical sixth-order forward algorithm. However, Omelyan, Mryglod and Folk^ have derived a number sixth-order 
gradient based, but non-forward, algorithms of the form 

J-jfi) ( £ ) = . . . e e{v a V+e 2 u„U) e et 1 T e e(v 1 V+e 2 u 1 U) e et 2 T e s(v 2 V+e 2 u 2 U) e et 3 T _ ^ 

Since the factorization is left-right symmetric, we only list the operators from the center to the right. These sixth-order 
algorithms all require a mixture of five force or effective force evaluations. Figure 7 shows the convergence of four of 
their position-type, six-order algorithms OMF40, OMF41, OMF43 and OMF45 when solving the Coulomb potential. 
None exhibited sixth order convergence. The best is OMF45, which converges with power 4.5, same as the optimized, 
supposedly sixth-order algorithm 4C opt with a — 0.41. Why these sixth-order algorithms are so down-graded in the 
Coulomb case is not understood. 

In Figure 8, we compare all algorithms on an equal effort basis as discussed earlier. In this case, the optimized 
fourth-order 4C algorithm has the smallest error, even when compared with OMF's sixth-order algorithms. 

In Figure 9, we show the convergence of these four sixth-order algorithms in solving for the ground state energy 
of the spiked harmonic oscillator with M = 6 and A = 0.001. All OMF algorithms can now be well-fitted with sixth 
order power-laws as indicated by solid lines. In the case of OMF40 and OMF41, the "glitch" in the convergence 
curve near e = 0.011 is real. The convergence curve for these two algorithms contains a singular term of the form 
w l/(e— 0.011), which blows up near e 0.011. Why only algorithms OMF40 and OMF41 exhibit such a singular 
behavior is also not understood. 

In Figure 10, we compare these gradient algorithms in an equal effort basis. The convergence range of sixth-order 
algorithms are not longer than those of fourth-order algorithms. For e > 0.002, the optimized fourth-order algorithm 
4C with a = 0.22 has smaller errors than all the sixth-order algorithms. However, for very high accuracy, sixth order 
algorithms are better when e is very small. 
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VII. CONCLUSIONS 



In this work, by regarding the radial Schrodingcr equation as a classical time-dependent force problem, we have 
shown that the entire literature of symplcctic integrators can be used to find its solution. Among symplectic inte- 
grators, factorized algorithms are favored because Suzuki's rule can be applied easily to solve the time-dependent 
force problem. Among factorized algorithms, gradient or forward algorithms are particularly suited because they 
take advantage of the harmonic character of the Schrodinger equation. We demonstrated the unique effectiveness 
of fourth-order gradient symplectic algorithms in solving the radial Schrodinger equation via Killingbeck's backward 
iteration. Even for very singular potentials, these algorithms are highly effective in computing the eigenvalue-function 
pair. There is also no difficulty in obtaining excited states. These gradient algorithms can form the core basis for 
solving non-linear Schrodingcr equations such the Hartree-Fock and the Kohn-Sham equation. However, due to the 
unique identification of the one-dimensional spatial coordinate as time, the current method does not appear to be 
generalizable to higher dimension for solving the general Schrodinger equation in two- or three-dimension. 

Among gradient algorithms, algorithm 4C with a tunable parameter a is most efficient in solving a variety of 
different potentials. Despite the fact that there are more complex fourth or sixth-order algorithms which use more 
effective force evaluations, none are really better than 4C. More force evaluation docs not necessarily enhance the 
efficiency of an algorithm, specially in solving the radial Schrodingcr equation. 

In solving the Coulomb potential, some gradient algorithms are down-graded to lower order while others are not. 
Even more surprising is the fact that none of the sixth-order algorithm exhibited sixth-order convergence in the range 
of e considered. These findings are not understood and should be studied further. 

By regarding the radial Schrodinger equation as a classical dynamical problem, one can now use the same set of 
symplectic algorithms for solving both classical and quantum mechanical problems. 
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TABLE I: Equal computational effort comparison of all fourth order algorithms. 
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SHO(M=4) 
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£ = 0.0002 


4B(e) 
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1.53438158545 


1.502005626 
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FIG. 1: The convergence of various fourth order algorithms in solving for the ground state energy of the Coulomb potential. 
Solid lines only connect data points to guide the eye. See text for identification of each algorithm. 




FIG. 2: Equal computational effort comparison of selected algorithms in solving for the Coulomb ground state energy. See text 
for discussion. 
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FIG. 3: The convergence of various fourth-order algorithms in solving for the ground state energy of the spiked harmonic 
oscillation 1)5. 2|l with M — 6 and A = 0.001. Same plotting symbols are used to designate the same algorithm compared in 
Figure 
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1.6399280 




FIG. 4: Equal effort comparison of various fourth-order gradient symplectic algorithms in solving for the ground state energy 
of the spiked harmonic oscillator of Figure El 
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FIG. 5: The convergence of various fourth-order algorithms in solving for the ground state energy of the spiked harmonic 
oscillation E2J witn M = 4 and A = 0.001. 




FIG. 6: The convergence of various fourth-order algorithms in solving for the ground state energy of the spiked harmonic 
oscillation 15,21 with M = 5/2 and A = 0.001. In this "supersingular" case, the power-law behavior is not observed within the 
range of e considered. 
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FIG. 7: The convergence of various fourth and sixth order gradient algorithms in solving for the ground state energy of the 
Coulomb potential. The solid lines here are fitted power-laws of power 3.5 (4B), 4 (4C, OMF40, OMF41, OMF43), and 4.5 
(4C with a = 0.41, OMF45). None is showing sixth order convergence. 




FIG. 8: Equal effort comparison in solving for the ground state energy of the Coulomb potential. The OMF algorithms are 
nominally sixth-order algorithms. However, their convergence is no better than that of algorithm 4C with a = 0.41. 




FIG. 9: The convergence of various fourth and sixth-order algorithms in solving for the ground state energy of the spiked 
harmonic oscillator (15.211 with M = 6 and A = 0.001. The solid lines are fitted power-laws of power 4 (4B, 4C), 5 (4C with 
a = 0.22), and 6 (all OMF algorithms). 
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FIG. 10: Equal effort comparison of various fourth and sixth-order gradient symplectic algorithms in solving for the ground 
state energy of the spiked harmonic oscillator of Figure [5] The optimized algorithm 4C with a = 0.22 has the smallest error 
for e > 0.002. 



